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Introduction 


The  purpose  of  this  investigation  is  to  find  theoretical  solutions 
to  ottering  of  electromagnetic  fields  by  burled  obstacles.  The  field* 
are  generated  by  wits  antennas  above  the  ground  and  the  burled  obataclee 
may  be  either  conductor  or  dielectric  bodies.  The  effort  of  thla  investi- 
gation is  divided  into  three  major  parts  as  follows: 

(I)  Wire  antennas  above  the  lossy  ground 

(II)  Scattering  by  a burled  wire  conductor 

(III)  2-dlaensional  scattering  by  dielectric  bodies. 

The  first  part  of  the  investigation  is  primarily  concerned  with  the  source 
problem  of  the  ccactering  phenomena.  It  enables  us  to  find  the  near  fields 
of  an  antenna  above  a lossy  ground  and  the  near  fields  penetration  into 
the  ground,  which  are  the  primary  fields  to  be  scattered  by  the  obstacles. 
In  order  to  study  that  problem  the  Somserfeld's  integrals  are  extensively 
studied.  The  integrals,  which  are  part  of  the  classical  solution  to  the 
Green's  function  of  Maxwell's  equations  Involving  lossy  half  space,  are 
useful  Co  our  problem  only  if  they  can  be  computed  economically.  Hence, 
the  first  part  of  th..s  f.  ivestlgation  is  the  extensive  search  for  techniques 
which  speed  up  the  convergence  of  Sommcrfeld's  integrals. 

The  second  part  of  the  investigation  involves  some  practical  applica- 
tion of  the  results  of  the  first  part  which  have  direct  bearing  on  the 
"detectability"  of  buried  obstacles.  It  may  answer  some  fundamental 
questions  relating  the  input  power  of  the  primary  source  to  the  scattered 
fields  from  a buried  obstacle  back  to  the  air. 

The  third  pj.rfc  of  the  invest igation  contains  the  preliminary  work  of 


scattering  by  burled  dielectric  obstacles.  Since  the  \inal  objective  of 
this  research  Is  the  detection  of*burled  dielectric  bodies,  the  technique 
of  computing  dielectric  scattering  is  of  primary  importance.  Unfortunately, 
there  has  been  no  such  work  reported  in  the  open  literetures  in  field 
theories.  Hence,  we  have  to  undertake  the  investigation  of  the  most  basic 
work  in  dielectric  scattering  problem  here.  Our  work  on  "Unl-moiaent 
method"^  has  advanced  the  computational  methods  to  Include  scattering 
by  dielectric  bodies.  Hie  application  of  the  "Unl-moaant  method"  to  2~ 
dimensional  scattering  problems  is  presented  in  this  report.  It  has 
demonstrated  that  the  "Uni -moment  method"  indeed  solves  many  dielectric 
scattering  problems  which  are  hitherto  unsolvable.  The  work  reported  in 
this  investigation  has  layed  the  ground  work  for  further  study  of  scattering 
by  burled  dielectric  bodies.  The  prospect  of  solving  a problem  which  haa 
been  thought  tn  be  untouchable,  is  better  than  ever. 

"This  rewearch  was  supported  by  the  Mine  Detection  Division, 

Countermine /Count or  Intrusion  Department,  U.S.  Array  Mobility  Equipment 
Research  and  Development  Center,  and  by  the  National  Science  Foundation." 


PART  OWE 

WIRE  ANTENNA  AliOVE  THE  LOSSY  GROUND 


PART  TWO 

SCATTERING  BY  A BURIED  WIRE  CONDUCTOR 
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I.  Fiald  In  fir  ala  and  Antenna  Intagra.1  Equation* 

Elementary  field  theory  stetee  that  the  electric  (laid  dua  to  an 


infinitesimal  currant  element  oriented  in  tha  a -direction  ia 

atU)  ■ *<'•)  1 >7  e !>I  c<^'>  - J““c<r'r’>*!  2r  • 

-Jkjr-r'l 

where  G(r,r')  * - — — (or  the  fraa  apace, 

I r-  r ( 


[2] 


d-1) 


and  r',  r are  the  position  vectors  of  the  source  and  field  respectively. 
Thus,  the  electric  fields  of  a linear  antenna  of  length  L ia  the  integral 
of  eq.  (1-1> 


E(r)  - ~ f G(r.r')  - ju>pG(r,r’)i]  I(r')d£’  . 

L J 


(1-2) 


Therefore,  when  the  current  1(a)  on  a linear  antenna  ia  known  it  ia  only 
necessary  to  carry  out  the  integral  of  aq.  (1-2)  to  find  the  electric 
fields  everywhere  in  apace.  The  fundamental  problem  ir.  antenna  theory 
is  to  fiud  the  current  due  to  a specific  excitation. 

Eq.  (1-2)  la  an  integral  representation  of  the  electric  field*  in 
terms  of  the  current.  If  we  let  r approach  to  the  conducting  surface  of 
the  sntenna,  where  the  tangential  electric  field  E (r)  Is  known,  i.a., 


E (r)  - - E*"C(r)  , 


(1-3) 


wher>.  E*nc(r)  is  the  incident  electric  field;  "Jq.  (1-2)  becomes 


/ lA-G(r.r’)  + k2G(r,7*)]  I(r’)d*’  - - j w c E;*‘v(x). 

L 3x*  * 


.inc, 


(1-4) 


Thv1!  above  equation  is  an  integro-dif ferential  equation  from  which  the 


Eq.  (1-4)  represents  en  sntenna  problem 


current  I(s')  mey  be  found. ^ 

Inc 

If  (s)  is  s localised  f<eld  such  as  that  produced  by  a generator,  and 

Inc 

it  represents  a scattering  probleai  it  E ^ (t)  is  an  illuminating  field 

such  es  a plane  wave. 

Thu  above  bricl  discussion  In  limited  to  linear  antenna  or  acattvrer 
in  open  uniform  space,  because  G(r,r'),  as  given,  is  the  solution  of  the 
scalar  wave  equation  in  an  infinite  uniform  space.  The  method  ie  readily 
extended  to  antennas  or  scatterars  in  or  above  a lossy  half-space  if  the 
Green's  function  is  replaced  by  the  solution  of  a acalar  wave  equation 
satisfying  the  boundary  conditions  of  the  lossy  hslf-apaco.  Such  Green's 
functions  are  known  a»  the  Sommerfeld's  integrals ^ represented  by 
*(r,r').  The  details  of  derivations  of  ths  Sommer fold's  integrals  may 
be  found  in  reference  (4)  and  given  in  ',/pendix  A,  so  we  shall  not  rederive 
the  integrals  here. 

In  principle  all  one  has  to  do  in  to  substitute  the  proper  Soamerfeld's 
formula  in  place  of  G(r,r')  in  (1-4)  in  order  to  find  the  current  on  a 
vertical  or  horison^.al  linear  antenna  above  a lossy  half-space.  Such 
applications  result  in  immediate  dlf faculties  because  Sommerfeld's  Integrals 
converge  very  slowly.  For  instance,  ir;  order  to  solve  a practical  antenna 
p-oblem,  the  Green's  function  G(r,r')  has  to  be  evaluated  for  different 
values  of  r and  r'  At  least  a few  hundred  times.  The  evaluation  of  any 
one  of  the  Soomier feld 'a  integrals  will  require  seconds  for  esch  pair  of 
values  of  r end  r',  il  we  use  a conventional  integration  routine.  Thus, 
the  application  of  the  integral  equation  is  economically  unfeasible,  unless 
we  can  increase  the  speed  of  computation  of  Sommerfeld's  integrals.  In  i u 
following  chapter,  we  shall  discuss  several  methods  of  approximation  of 
Sommerfeld's  integrals  for  fast  convergence. 


1 1 . Approximation  of  Sommerfald's  Integrals  for  Tea t Convergence. 

SoMMrfeld'a  Integrals,  the  Graar.'e  functlona  of  Maxwell's  equations 
eatlsfylng  the  boundary  conditions  of  lossy  haH-sparc,  were  first  derived 
by  A.  Sommer  1 eld  in  1909.  These  integrals  have  since  been  extensively 
studied  by  many  scientists  in  search  for  the  ground  waves.  In  the  applica- 
tion of  Soamerfeld's  integrals  to  study  the  antennas  over  loesy  ground,  we 
are  not  so  much  Interested  in  the  physics  of  the  problem  as  we  do  the 
mathematical  or  computational  aspect,  of  the  integrals,  which  as  we  have 
ehovn  in  Chapter  I is  useful  only  if  the  integrals  converge  rapidly.  Un- 
fortunately, the  integrals  arc  oscillatory  and  slowly  decaying  type 
which  arc  most  time  consiming  to  compute. 

To  concentrate  our  attention  only  to  the  study  of  the  convergence  of 
the  Soamerfeld's  Integrals,  we  have  left  the  derivation  of  the  integrals 
in  Appendix  A.  Only  the  derivation  of  the  Integrals  for  V'e  horlxontal 
dipola  antennas  are  given  in  the  Appendix,  because  the  formulas  of  a vertical 
antenna  are  essentially  contained  in  those  of  the  horlxontal  ones. 

The  Soamerfeld's  integrals  of  an  horixonal  dipole  antenna  are  gi\ on 
by  (A-ll)  - (A-16)  of  Appendix  A,  but  are  given  in  the  following  also  for 
convenience: 
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Vg*-vh  v 
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(2-4) 


where,  the  source  Is  assumed  to  be  in  the  air,  n and  tt  are  rc:  tively 

the  x and  2 components  of  the  Herts  potential  meas  *c  ' . and 

it  ..  and  x ore  the  sane  components  measured  .he 

Xr.  ZE 

complex  refractive  index  of  the  ground,  ano 


2 2 
A - k 


VE  " kE  • 


2 J j u2  e 

j and  k^  • — j—  are  the  wtve  numbers  in  air  and  earth  respectively  , 
c c 

The  constant  h is  the  height  of  the  source  location  and  the  variable  z is 
the  height  of  the  field  location  as  shown  in  Fig.  A-l. 

When  the  source  is  located  in  the  ground,  ue  shall  interchange  and 
71  R,  v and  Vg,  and  replace  k^,  ^ for  k,  1 in  (2-1)  and  (2-2),  and  similar 
modifications  for  (2-3)  and  (2-4). 

The  approximations  we  shall  use  take  advantage  of  the  decaying  property 

of  the  exponential  which  is  quite  different  for  (2-1)  and  (2-2),  and 
v(z+h)  Vg(z+h) 

whether  it  is  ev^z  or  e . Therefore,  it  is  necessary  that  the 

approximation  techniques  are  treated  in  3 different  cases,  namely,  (A) 
both  source  and  observation  locations  above  the  ground  (h  > 0,  z > 0), 

(8)  source  location  above  and  observation  location  under  the  ground 
(h  -*  0,  z c 0),  and  (C)  both  source  and  observation  locations  under  the 
ground  (li  0,  z ^ 0) . 
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The  approximation  technique*  of  th>  case  are  numerous.  They  are 
reported,  for  instance,  in  references  [6]  and  [7].  They  have  asdumed, 
hovever,  the  refractive  index  of  the  half-space  to  be  large  and  the 
observation  point  to  be  far  away  from  the  Interface.  In  this  section 
wo  sh.iJ  I use  an  .ipproxlmalc  boundary-condition  technique  to  arrive  to 
an  approximation  ot  Sommuriuld's  integral,  which  is  ajso  valid  for  source 
or  field  points  relatively  close  to  the  ground. 


roxixate  Boundary  Condition 


Consider  the  Hertz  potential  1rx(°r  *z)  over  a small  portion  of  air- 
earth  interface.  The  dependence  of  all  potential  components  on  the  radial 
direction  may  be  described  by 


ir  - w exp(-J  (k^x  + kyy)),  z > 0, 


(2-5) 


where  w is  the  so  called  "attenuation  factor"  and  is  a function  of  position. 
It  is  assumed  that  the  attenuation  factor  w does  not  vary  too  much  over 
the  range  of  a wave  length  Aq  in  the  air,  or  quantitatively. 


(2-6) 


X 

A . 

_Jo_  aw  < , -l . 

o 3w 

2wu)  3x 

2iru>  3y 

I 8 I 

This  assumption  is  claimed  valid  in  most  of  the  practical  situations. 
With  equation  (2-5)  and  (2-6) , the  solution  of  the  wave  equation 

2 2 

a £ 


V + kgiTg  - 0 for  z 5 


(2-7) 


is  obtained: 


where  n is  the  Hertz  potential  for  the  source  in  the  air  and  observer  in 
E 

2 2 2 

the  ground,  and  k ■ k + k . If  equation  (2-3)  is  substituted  into  the 
a x y 

boundary  cordition  for  it  at  z "0,  (eq.  (A-7)  and  (A-8)),  we  obtain  a new 
set  of  boundary  conditions  valid  for  the  solution  in  the  region  z > 0 


ii  at 
x 


z 


0 : 


(2-9) 


3tt 

* 

3x 


at  z 


0 , 


(2-10) 


where 


o 2 

- cos  ip  , 


(2-11) 


and  e°  * c2/tj.  (2-12) 

The  ip  in  eq.  (2-10)  is  defined  according  to  the  propagation  constant  and 
is  given  by  cos  ip  ■ k^/k.  Since  k is  real,  ip  is  real  and  hereon  will  be 
considered  equal  to  the  angle  that  a wave  is  incident  on  the  zeroth 
Fresnel's  zone.  The  above  assumption  is  justified  by  the  following  con- 
sideration. From  the  theory  of  wave  propagron  between  transmitter  and 
receiver  above  a ground  plane,  it  is  found  that  the  wave  which  arrives 
at  the  receiver  from  the  contribution  of  reflected  wave  is  essentially 
from  one  direction.  That  is,  there  is  only  a certain  limited  area  on 
the  ground  plane  from  which  a reflected  wave  may  reach  the  receiver. 

Usually,  for  practical  applications,  only  a few  Fresnel's  zones  are 
needed  to  be  a good  approximation.  The  only  case  when  the  wave  zone  covers 
a large  area  is  when  both  the  source  and  observer  are  very  close  to  the 
ground.  But  for  this  case  cos  ip  is  almost  equal  to  1.  Therefore  our 
assumption  of  optical  reflection  angle  is  practically  valid. 
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The  approximate  boundary  conditions  of  eq.  (2-9)  and  (2-10)  are 
similar  to  Che  impedance  condition  except  that  the  "effective  refractive 
indecea"  are  used  here.  With  these  boundary  conditions,  we  obtain  a set 
of  solutions  to  the  wave  equation  if  we  follow  the  similar  procedure  as 
that  in  Sommerfeld's  treatment ^ ^ . The  solutions  are  given  by 


-JkR  -jkR'  -jkR' 

e e ..  {2  e 
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R' 


R' 


O 

- 2jkVe^  f 


J0(Ar)e 


-v(z+h) 


jkv  t 


dA 


o v 


OD 


(Ar)e 


-v(,+h) 


A dA  , 


1 

(2-13) 

(2-14) 


e v + jityio 

where  cos  4 - x/r,  R - >4r2  + (z-h)2,  and  R*  « J r 2 *•  (z+h)^ . The  Integrals 
in  the  above  formulas  can  be  transformed  to  other  integral  forms  which  do 


not  involve  the  Bessel's  function.  The  detail  of  the  transformation  is 
given  in  reference  [9].  In  Appendix  B,  the  integral  in  Eq.  (2-13)  is 
shown  to  be  equivalent  to  the  following  integral. 


z 

I exp(-jk[R£  + (£-z)])d£/R^  # 

> 


The  Hertz  potentials  are  then  given  by 


,-jkR  -jkR'  -jkR' 


R R' 


{2  ; 2jk  y/c? 
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rJ°°  _ 


jk[[R’  + /TjU-z)] 
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(2-15) 


*z  “ k2  dx  **A  + V ’ 


(2-16) 


where  1^  and  1^  are  given  by 
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and  R' 

V» 


+ <h+02. 


For  the  convenience  of  comparison  let 


The  equation  U is  corresponding  to  the  integral  term  of  Sooner f eld's 
formula  in  equation  (A-ll)  i.e. 


U 
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2 f Jo(Ar)e-',(t+h) 
JQ 


X 

v + vE 


dA. 


The  equation  of  V can  be  shovn  to  be  corresponding  to  the  Sommerf eld's 
formula  for  the  vertical  dipole  and  is  given  by  [3] 


-jtt’ 
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Numerical  Results 


The  integrals  in  U and  V of  eqs.  (2-19)  and  2-20)  are  computed  by 
using  the  Simpson  algorithm  in  a complex  plane.  The  choice  of  integral 
path  is  rather  arbitrary.  We  find  that  a straight  line  path  extending 


from  £ • x io  £ ■ x - J®  is  good  enough. 

The  values  of  U and  V are  given  in  Pigs.  (2-1)  to  (2-6)  for  lossy 
tangent  being  .01  and  1 and  different  locations  of  observer  and  source. 
Values  of  and  Vg  in  Eqs.  (2-21)  and  (2-22)  are  also  computed  and  shown 
in  these  figures  by  dots.  It  is  clear  t l*.-t  the  approx  ini  it  ion  given  by  U 
and  V is  accurate.  The  computational  time  we  consume  is  of  the  order  of 
100  times  less  than  that  using  U and  V . Together  in  these  figures  the 

8 S 

familiar  approximations  using  the  plane  wave  reflection  coefficients  are 
also  presented.  As  expected,  the  plane  wave  approximation  is  not  a good 
approximation  where  value  of  (h  + *)  is  getting  smaller. 

Taylor's  Expansion  of 

The  approximation  of  (A-ll)  and  (A-15)  by  (2-13)  and  (2-14) , while 
physically  justified  in  previous  paragraphs,  is  mathematically  amount  to 
replacing  v£  by  jk  /t'“,  which  makes  the  transformation  to  the  integrals 
of  (2-15)  and  (2-16)  possible.  Approximating  v^  - \A\2-k^  by  a 

constant,  in  an  Integrand  must  also  be  justified  mathematically.  Indeed, 
it  can  be  and  we  shall  show  £.n  the  following. 


Let 


^(X2-k2)  - (4-k2) 


J~1 

» vv  - 


2 T 


(2-23) 


The  Taylor's  expansion  of  (2-23)  at  v * 0 is 


; j /k2-k2  (1  + ( *k?-k2)v2  + higher  order  terms) 


i*.f\  (1  [1+  k/^u  -isls!*  )v2+  ••■]  • (2-24) 
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DISTANCE  IN  ME'kERS) 


the  U integral  vs.  horizontal  range,  h > 
Where  R:  plane  wave  approximation,  AST: 


boundary  value  method. 


— 

0 


0.1  02  0.3  0.4  05  0.6  0.7  0.8 

r (HORIZONTAL  DISTANCE  IN  m) 


Fig.  2-3.  Value  of  the  U integral  ve.  horizontal  range,  h > 0, 
* > 0 , (REF  la  equivalent  to  R,  and  R and  ABT  are 
defined  in  Fig.  2-1.) 
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Fig.  2-4.  Value  of  V integral  va.  horiiontal  range  h > 0,  x > 0. 
(R£F , ABT  are  defined  in  Fig.  2-3.) 
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It  is  noticed  that  v„  ? Jk  /]*  ia  Just  the  first  term  of  the  expansion. 

The  expansion  of  (2-24)  is  valid  for  |v  sin  ♦ ( • 1,  which  is  violated  when 
| A | > > { k | . However-  when  that  happens,  the  exponential  tera  e 
will  be  sufficiently  snail  to  suppress  the  integrand.  That  is  to  say  the 
value  (h*s)  has  to  be  sufficiently  large  to  aake  Che  suppression  possible. 
For  frequency  at  1 ('.Hr,  we  find  that  the  value  of  (x+h)  can  be  as  sauill 


as  0.05  outer. 


B.  (h  > 0,  x - C)(10) 

The  SoMierfeld's  formulas  are  given  by  eqs.  (2-2)  and  (2-4).  There 

is  v appearing  in  the  exponential  tera,  which  changer  the  behavior  of 

the  above  integrands  drastically.  Now  the  exponential  tana  starts  to 

decay  for  auch  larger  value  of  X than  that  in  case  A.  Therefore  a auch 

better  approximation  technique  is  required. 

The  discussion  of  Taylor's  expansion  of  v>  in  previous  section  leads 

ua  to  believe  that  aore  terms  in  the  expansion  should  be  adopted  for  a 

-vEh 

better  approxiaation  of  e Therefore  the  following  expansion  is 

considered. 


V 


U + 


2 

v 

2! 


(2-25) 


The  above  foraula  is  valid  only  when  v is  saall.  When  X Increases  toward 

Ke(k„),  the  right  hand  wide  of  (2-25)  grows  quadraticelly , while  the  left 
c 

hand  side  is  actually  a decaying  function.  Furthermore,  the  quadratic 
giowth  is  rather  rapidly  before  the  e is  able  to  suppress  the  integrand, 
contributing  a substantial  amount  of  error  to  the  integration. 

Aiming  at  eliminating  the  quadratic  growth  of  the  expansion,  we  then 


com  up  with  the  following  approximat ion 

v -V'W) 


2 *v  (x  v) 

( I ♦ * -,)  ♦ »■  " - ( I ♦ t v + — ) (2-26) 

/ 2 2 " ° 

t 


X V 

Where  in  (2-26)  we  have  purpoaely  added  t equation  (2-25)  a term  e ° and 
aubtracted  its  second  order  Taylor's  expansion.  In  order  to  eliminate  the 

second  order  term,  z is  so  chosen  such  that 

o 


(i  v)2 
o 

2 


«xp(-J 


*e 


«)/[«v*/2j 


'/kfk2)] 


(2-27) 


The  expansion  of  (2-26)  ia  now  only  of  linear  growth.  From  Table  (2-1) 
we  can  see  that  considerable  improvement  has  been  achieved  by  using  this 
technique . 

One  of  the  limitations  makes  the  approximation  so  difficult  is  that 

functions  that  are  used  to  approximate  it  must  be  functions  of  the  form 
XV 

e This  limited  forms  of  functions  ensure  the  further  reduction  to  a 

V 

fast  convergent  integral  possible.  From  the  plotting  of  e with  respect 

V 

to  A,  we  find  that  the  behavior  of  e is  rather  smooth  when  A is  small. 

2 

This  is  because  when  A is  small  the  effect  of  A is  neglegible  compared 

2 VEZ  -vh 

with  k . As  A is  getting  larger,  e becomes  oscillatory;  but  then  e 

is  small  enough  to  suppress  the  contribution  of  this  portion  to  the  whole 

integral.  Hence  a better  approximation  may  be  obtained  if  we  have  some 

V 

way  to  approximate  the  function  e while  A is  not  too  large.  This  brings 
the  idea  of  curve-fitting  method. 

With  the  experience  oi  pievious  approximation  methods  which  led  to 
eq.  (2-26)  we  conclude  that  a reasonable  formula  for  the  curve-fitting  is 
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VALUE  OF  A 


L ftlltl 


.u*  ,-.¥i 


Fig.  2-9.  Comparison  of  curve-fitting  approximation  and  exact  formula 
f = 109Hz,  f = 15-j. 1592 , z = -.45m. 


A is  defined  in  Fig.  2-7. 


In  (2-28),  zq  is  determined  the  same  way  as  that  in  (2-27),  and  the  co- 
efficients are  found  by  enforcing  Fr  - F(A)  at  four  points,  namely,  A ■ Aq, 


*2*  *3*  The  choice  of  depends  upon  k,  k^  and  z.  Usually,  two 

A^s  are  chosen  in  A < k,  one  at  A - k and  the  other  at  around  (l/2)k£. 
Note  that  eq.  (2-28)  differs  from  eq.  (2-26)  by  having  more  freedom  on 
the  weight  of  each  term.  The  zeroth  order  term  of  Taylor’s  expansion  in 
eq.  (2-26)  is  included  in  a^  of  eq.  (2-28). 

For  the  comparison  of  Fn  and  F(A),  values  of  the  integrand  of  eq. 
(2-2)  are  computed  for  e » 15  - J*1592,  f » 1GHz,  and  ore  shown  in  figs. 
(2-7)  to  (2-10).  It  is  seen  that  Fn  stays  to  be  a good  approximation  up 
to  | z | > .45  meter.  From  Fig.  (2-10),  with  z ■ .6  meter,  or  equivalently, 
two  free  space  wave  lengths,  the  approximate  curve  fails  to  respond  to 
the  oscillatory  behavior  of  F(A)  when  A > 30.  This  failure  when  z is 
so  deep  under  the  ground  plane  is  not  a serious  matter,  because,  the  value 
of  potnetial  is  so  small  then,  that  any  approximation,  including  plane 
wave  approximation,  tends  to  converge  to  the  correct  value. 

With  the  approximation  of  eq.  (2-28),  we  are  able  to  transform  eq. 
(2-2)  into 


"jcE  ■ 4 ‘V, 


+ Vo  + Vi 


V2> 
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where 
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J (Ar)e 
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-v(h-z) 
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Jk 


AdA 
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(2-29) 


(2-30) 
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Wherein  the  denominators  of  the  above  Integrals  we  have  replaced  vg  by 

jk  ^ as  we  have  done  in  previous  sections.  By  a similar  procedure  as 

that  in  section  A where  a > 0,  I is  reduced  to: 

o 


(2-34)  ] 


where  R"  - vr*  + h2  and  R’.’  - >42  + £2.  From  eqs.  (2-32)  and  (2-33),  j 

* I 

1^ , are  easily  obtained;  their  final  forms  still  contain  the  integral  j 

in  (2-34).  1 in  eq.  (2-30)  is  obtained  by  replacing  h in  eq.  (2-34)  | 

1 

by  (h-s  ).  i 

As  a comparison  of  accuracy  by  using  (2-25),  (2-26)  and  (2-28)  for  i 

the  approximation,  values  of  are  tabulated  in  Table  (2-1).  Figs. 


(2-11)  and  (2-12)  show  the  value  of  ir  in  eq.  (2-29)  as  a function  of 

xE  „ t 

Table  2-1  Comparison  of  Different  Approximation  by  Using  £qs.  (2-25) 
(2-26),  (2-27)  and  Exact  Function. 


'•.mi*-,.  > 


Nv  "xE 

z(«o\ 

exact 

curve  fitting 

eq.  (2-24) 

eq.  (2-23) 

-.1 

-.0977  4-  j . 1130 

-.0998  + j .1142 

-.1061  4-  J.1110 

-.1239  4-  j.130 

-.2 

.0694  j .0376 

.06912  + j .0363 

.0639  4-  J .0433 

.0895  + j .082 

z (f ixed  r)  and  function  of  r(fixedz)  at  iCIfc  and  » *=  15-J  1.592.  Exact 
formula  of  eq.  (2-2)  is  also  evaluated.  It  is  seen  from  these  figures 
that  the  curve-fitting  method  yields  a very  good  result.  It  is  well  known 
u*.at  the  wave  propagating  inside  a lossy  ground  from  a source  in  the  air 
behaves  like  a plane  wave.  Hence  we  also  compare  our  curve-fitting  method 

with  plane  wave  approximation  and  ahow  the  result  in  fig.  (2-13).  Only 
the  m itudes  are  shown,  because  the  phases  agree  very  well  everywhere. 

It  is  observed  that  while  the  phase  behaves  well  the  plane  wave  approxi- 
mai  may  result  an  error  up  to  40*  in  magnitude.  Earlier,  we  have 
obser  from  fig.  (2-10)  that  when  z gets  deeper  in  the  ground  the  accuracy 
will  deteriorate.  The  degree  of  inaccurac., , however,  is  not  significant 
judging  i.om  fig.  (2-13).  Besides,  the  value  of  is  sufficiently 

small  wt,_..  Inaccuracy  is  observed.  In  that  case  the  information  of  phase 
is  more  important  than  that  of  magnitude. 

The  ( proximation  of  eq.  (2-26)  can  also  be  applied  to  eq.  (2-4)  to 
obtain  an  approximation  to  w . Then  eq.  (2-4)  becomes 


’zE  * 72  h [Vz  * Vo  + Vl  + V2J 

KE 


(2-35) 


Where 


-v(h-z  ) v - jk  %/c^ 

J (Ar)e  0 AdA  , 

° v 4-  jk/^ 


(2-36) 
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(2-37) 
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Jo(Ar) 


T=rc  AdA 

cv  + Jk  ve° 


.»2, 


)2h2 


(2-38) 


(2-39) 


Comparing  (2-37)  with  (2-14)  we  find  that  an  expression  of  the  fora  (2-16) 
for  Iq  in  (2-37)  is  immediately  obvious. 

C.  (h  < 0.  s < 0) 

The  Sommerfald' a formulas  are  now  given  by 


* 
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- J4Rt 

R* 
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J0(Ar) 


-v£(*+h) 
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v£  ♦ v 


dA 


(2-40) 
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^2 
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T cos  4 
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Jj(Ar)e 


-vE(s+h) 


- v 


Cv  4-  V, 


A2dA 


(2-41) 


The  integrals  in  (2-40)  and  (2-41)  are  siailar  to  those  in  eq.  (A-ll) 
and  (A-1S)  for  the  case  A except  now  we  have  v_  instead  of  v.  Because  of 

b 

this  change  we  have  to  approximate  v in  the  denoainators  of  the  above 

integrands.  Approximating  v by  a constant  is  obviously  not  satisfactory, 

since  the  value  of  v varies  substantially  in  the  range  0 < A < Re(kE>.  It 

is  obvious  we  have  to  look  for  an  approximation  of  v by  some  function  of 

A.  In  search  for  this  suitable  function  of  A we  also  have  to  bear  in  alnd 

that  the  approxiaate  function  is  eventually  to  be  coupled  with  v and 

£ 


made  the  transformation  of  integral  into  a faat  convergent  fora  possible. 

The  experience  of  previoua  sections  leads  us  to  try  the  approximation 
1 n ci 

of  — by  £ . Therefore  we  look  into  the  possible  approxiaation 

V VE ' i VE  di 

of  v by  a polynomial  of  v_.  The  tera  l/(v+  v,. ) then  can  be  reduced  to  a sub 

c*  b 

of  terms  of  the  fora  1 / ( v>.  + constant)  by  partial  fraction.  Theoretically, 
polynoaial  approxiaation  of  a function  is  feasible  in  an  analytical  region. 
But  we  are  restricted  by  the  requirement  of  coaputing  efficiency  such  that 
only  a few  terms  (hopefully  less  than  three)  of  l/(v  + constant)  can  be 

b 

to  approxiaate  l/(v  + v„) . Our  experience  concludes,  however,  polynoaial 
approxiaation  of  v by  v_  is  not  satisfactory. 

XL 

The  behavior  of  rational  function  composed  of  two  polynomials  Is  more 
versatile  than  that  of  a single  polynomial.  While  keeping  the  order  of 
each  polynoaial  down,  rational  function  has  aore  degree  of  freedom.  In 
order  to  appreciate  this,  let  us  consider  the  following  approxiaation: 


A * ve  + *< 


(2-42) 


Therefore , 

__1 

v + vr 


Ve  + Ve  + *6 


(1  + a4)vj!  + <ax  + m5)s>l  + <a2  + a^  + >3 


(2-43) 


It  is  clear  now  that  we  can  then  separate  the  right  hand  side  of  (2-43) 
into  three  partial  fractions.  There  are  six  adjustable  constants  to  be 
determined  by  enforcing  the  equality  of  both  sides  of  eq.  (2-42).  This 
should  coapare  to  the  single  polynomial  approxiaation,  where  only  three 
adjustable  constants  are  possible  to  obtain  a sub  of  three  partial  frac- 
tions. 
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Equation  (2-42)  la  only  vclld  for  tha  ration  where  X >_  k,  ainca  v 

haa  a branch  point  at  X • k.  To  anaura  tha  continuity  of  v it  • 0 and 

correct  functional  bahavlor  whan  A extends  to  infinite,  wa  naad  to  fix  two 

constant!.  Nota  that  v ■ 0,  whan  A » k , and  *q.  (2-42)  becomes 

u E 


/*T-  ■ vu6 


. * . a^  ■ a^/  v kg  - k * A 


(2-44) 


Uhan  A -*■  •,  v • Vg,from  aq.  (2-42) 


A ^ - “4 


So  that  a^  should  ba  aqual  to  1,  and  aq.  (2-42)  becoaMs 


VE  ♦ Ve  * *2VE  + *3 


VE  + *5VE  + *6 


In  ordar  to  furnish  tha  approximation  for  eq.  (2-40)  lat 


(2-45) 


VE Ve  + Ve  + Ve 

y + VE  2Vg  + (ax  + *5) vg  + (*2  + *6)VE  + *3 


(2-46) 


The  right  hand  aide  and  left  hand  side  of  eq.  (2-46)  as  a function  of  A 
are  compared  and  shown  in  fig.  (2-14).  It  is  seen  that  a very  good 
approximation  beyond  A > k is  achieved.  For  that  particular  case, 

► » 15-J1.59  and  k - 20.958,  equality  of  eq.  (2-46)  is  enforced  at 
A - 25,  50t  100  and  a^  * 0.  Eq.  (2-46)  can  ba  decomposed  into  tha 
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+ s. 


(2-48) 


The  first  term  in  (2-48),  which  is  the  contribution  from  the  first  tern 
on  the  right  side  of  (2-47),  is  obtained  the  same  way  as  that  in  previous 
sections.  U^'a  are  given  by 


tt 

-■N 


— I JQ(*r)e 


-v£(z+h) 


AdA 


(vE  + VvE 


(2-49) 


which  is  the  similar  integral  to  the  integral  term  of  (2-13).  Hence 
can  be  easily  reduced  to  the  following  form: 
z 


„.!*/•  .-JVi  e-b.«-M>  a 

ui  2 / e e r; 

-j-  6 


(2-50) 


Again,  is  defined  the  same  as  that  in  section  A.  S in  eq.  (2-48)  has 
to  be  computed  by  brute  force  method. 
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(2-51) 
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f(v  ) is  given  by  eq.  (2-47)  and  the  integral  licit  is  X ; k. 

Ci  O 

Values  of  U^are  computed  according  to  (2-48),  and  are  shown  in  figs. 

(2-15),  (2-16).  Comparisons  are  made  between  results  of  eq.  (2-48)  and 

(2-40)  and  the  agreement  is  almost  perfect.  Computing  time  took  for  each 

point  of  calculation  is  less  than  .5  sec.  for  (2-48)  comparing  to  30  secs. 

for  (2-40)  at  a CDC  6400  computer. 

For  the  approximation  of  eq.  (2-41),  it  is  enough  to  approximate 
v£  - v 

in  terms  of  a rational  function  of  v„.  This  is  done  much  in  the 

cv  + vE  E 

same  fashion  as  that  we  have  just  presented.  With  eq.  (2-45), 
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(s  4 1)  j e + a5)  2 ^d2  £ + a6^  a3  £ 

E e 4-  1 VE  e + 1 E e 4-  1 


» (2-52) 


tne  right  hand  side  of  (2-52)  can  be  decomposed  into  partial  fractions: 
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V„  4-  b, 


4-  br 
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(2-53) 


The  comparison  of  both  sides  of  eq.  (2-53)  as  functions  of  X is  shown  in 
Fig.  (2-17).  All  what  follows  will  be  similar  to  the  procedures  we  have 


just  presented. 


'-'f* 
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Fig.  2-15.  Value  of  U£  (1/meter)  obtained  from  eq.  (2-48)  vs. 
horizontal  range  h < 0,  z < 0. 


Fig.  2-16.  Value  of  U (1/meter)  obtained  from  eq.  (2-48)  va.  vertical 


III.  Horisonta  1 Dipole  Antennas  Ovsr  Lossy  Half-Space 
Integral  Equation  Formulation 

a 

A a aentioned  in  Chapter  I that  the  antenna  problem  involving  loaey 
half-space  is  similar  to  that  in  free  space  with  the  free  space  Green's 
function  replaced  by  the  Sommerfeld's  integral.  In  this  chapter  we  £h«ll 
discuss  the  antenna  equation  in  greater  detail  and  the  Integral  equation 
formulation  will  be  described.  It  is  noticed  that  the  formulation  given 
in  Chapter  I is  an  lntegro-dlf ferential  equation  which  is  less  stable 
than  an  integral  equation. 

The  electric  fields  of  an  infinitesimal  current  can  be  expressed  in 
terms  of  the  vector  Herts  potential 

E - (k2  + V V •)*  , (3-1) 

where  satisfies  the  following  equation: 

(V2  + k2)+  - - 6(x-x')i  . (3-2) 

Applying  the  boundary  conditions  of  a lossy  half-space,  for  the  fields 
above  the  ground  (z  > 0),  41  happens  to  be. 

ij>  ■ t — (j»  x + n z)  , for  a x-oriented  dipole,  (3-3) 

HIT  X Z 


where  n and  « are  the  Sommerfeld's  integrals  given  by  (2-1)  and  (2-3) 

X z 

respectively. 

For  a horizontal  linear  antenna  oriented  in  the  x-direction,  we  have 
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4xjwr 


/ V‘'> 


(ff  x + n s)dv' , 
X z 


(3-4) 
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where  the  Integration  is  over  the  apace  occupied  by  the  antenna.  Substi- 
tuting (3-4)  into  (3**1)  we  get 


E(r) 


4nju»c 


f 


J (*') 

X 


(k  + VV 


) (n  x + 7t  z)dv' , 
x z 


(3-5) 


where  E is  now  the  electric  field  radiated  by  the  current  oi  the  antenna. 
Let  r,  the  field  point,  approach  the  surface  of  the  antenna,  the  total 
tangential  coaponent  of  the  electric  field  must  vanish,  thus 


E (r  ) + E (r  ) - 0, 
x'  a ox  8 * 


(3-6) 


where  Ex(rg)  i®  the  x-component  of  E(rg),  and  ^^(r^  i®  the  x-couponent 
of  the  electric  field  produced  by  a generator.  rg  is  a point  on  the 
antenna  surface. 

Since  the  x-component  of  V V • rr  is 


a2»  92. 

(vv  . .)  . i •<-?*+  ). 


3x 


(3-7) 


we  have,  after  substituting  (3-7)  and  (3-6)  in  (3-5), 


jw  Eox(x)  " h JIx(x')I(  ^7  + k2)\(x|x')  + 71  z (x | x ’ ) ]dx'  (3-8) 


where  a thin  linear  x-oriented  antenna  is  assumed,  and  I (x)  « 2naJ  (x) 

X X 

is  now  the  current,  and  a is  the  radius  of  the  thin  wire,  x and  x’  are 
points  on  the  antenna  and  L is  the  length  of  the  antei  ,ia.  Eq.  (3-8)  is 
the  in tegro- differential  equation  similar  to  Eq.  (1-4)  in  free  space.  In 
the  following  we  shall  convert  (3-8)  into  an  integral  equation,  which  is 
more  stable  for  numerical  computation. 
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It  is  well  known  that  the  differential  equation, 


( + k2)*(x)  - g(x) 

dx 


(3-9) 


has  the  aost  general  solution, 


* 5 i 


$(x)  * A sin  kx  + B cos  kx  + r I g(0  sin  k(x-Od£  • 


(3-10) 


It  is  obvious  from  (3-9)  and  (3-10)  that 

(~+k*)  ie  | 8(0  Sin  k(x-OdU  - g(x). 


(3-11) 


Now,  let's  multiply  both  sides  of  (3-8)  by  sin  k(x-0  and  integrate,  we 


get, 
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( -la,  r £ «)  sin  k(x-Od(  - — fl  (x')  (H  ( + k2)s  (t|x’)»in  k(x-Od{)dx’ 
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♦ h f V’*  <c  i k 1 ti  «>-OdC)  dx'.  (3-12) 
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u 0 


Using  (3-11),  we  find  the  first  term  of  the  right  hand  side  of  (3-12) 
is  simply 
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i / ‘i? 


+ k2)  TTx(5|x')sin  k(x-OdOdx' 


fcj 


Ix(x’)  nx^xl  x*)dx’  ” (A  sin  kx  + B cos  kx) . 


(313) 


Eq.  (3-12)  now  becomes 
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- A sin  kx  + B cos  kx 


- | EoxU)8in  k(x-OdC 


(3-14) 


The  second  tern  on  the  left  hand  side  of  (3-14)  contains  no  singularity, 

sc  in  confutation  we  may  differentiate  before  perform  the  integration. 

In  actual  computation  of  (3-14)  we  shall  use  the  fast  converging 

form  for  w and  ir  as  given  in  (2-11),  (2-1?)  and  (2-13).  It  should  also 
X z 

be  noticed  that,  although  (3-14)  involves  a double  integral  of  the  fora. 


P - 


/x 
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(3-15) 


we  can  decompose  (3-15)  in  the  following  terms. 


( lx(x’>  ' i (sin  kx 


a. 

I ac  ^ az 

X 
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"x(£lx')>  cos  kc  d£ 


- cos  kx  § — ( — *x(£|x'))  sin  k£  d£]  dx' } (3-16) 


In  (3-16)  the  integrands  are  no  longer  functions  of  x,  and  the  integration 
witli  respect  to  f.  to  diffe  values  of  x may  be  made  cummulatively  and 
avoid  redundant  computations. 


Numerical  Results 

Once  the  formulation  o ..  antenna  equation  is  properly  done,  the 
numerical  method  of  solution  is  rather  straight  forward.  Numerous  references 
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•re  available  on  numerical  technique*  [3],  [11].  In  hare  we  prcaant  a cm* 
typical  raaulta  of  computation.  In  Pig  (3-1)  the  Input  admittance*  of  a 
Horlsontal  dipole  antenna  na  a function  of  dlatance  above  the  lossy  ground 
•re  shown.  The  antenna  parameters  are  shown  In  the  figure.  It  Is  noticed 
that  the  proxlnlty  of  th'  lossy  ground  does  have  a visible  effect  on  the 
input  properties  of  the  antenna  and  thus  its  performance. 

The  prop, ram  we  have  written  is  also  capable  of  computing  the  horl- 
sontal antenna  arrays  connected  by  transmission  lines  such  as  shown  In  Figs. 
(3-2)  and  (3-3).  The  computed  results  we  have  presented  here  are  original. 
While  there  is  no  comparable  theoretical  or  experimental  results  available 
to  check  our  computation,  the  computations  of  a few  extreme  cases,  do 
indicate  that  the  program  is  free  of  bugs  and  gives  accurate  results. 

These  results  will  be  used  to  compute  the  incident  fields  for  the  analysis 
of  scattering  by  buried  obstacles  in  the  next  chapter. 
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ig.  3-1. 


Input  admittance  of  a horizontal  antenna  over 


a conducting  ground. 
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IV. 


Incident  .t ml  Scattering  Fields  for  l Ik*  Problem  of  Kurjed  Wire 


In  die  previous  chapter,  we  li.ive  discussed  the  method  ol  solving 
an  antenna  problem.  The  current  distribution  thus  calculated  will  be 
applied  hare  to  find  the  electric  field  everywhere.  The  field  equation 
(1-2)  will  be  restated  in  the  following  manner  for  an  x-oriented  current 
source: 


E1J(r) 


1 fV*"  2 

- J=^  J - V"  (kj 


+ U 


)(trlJx  + iriJ£)dl’ 

X z 


1,  2;  j - .1,  2 


(4-1) 


Where  doubla  indecea  are  adopted;  i and  j represent  the  locations  of 
source  and  field  point  respectively.  When  index  is  1 the  point  is  in 
the  air  and  2 in  the  ground.  Thus  » u *^uo*  kj  • 4 kj , and  Cj  « c , 

c2  * 1/c  • The  operator  V V • in  eq.  (4-1)  is  given  by 


V[V*(i»  x + x £)] 
x z 


2 2 2 

^ — [n  + $ ]x  + [ ■ - • ] [n  +■  + ]y  + [ — -r—  ) [n  + ♦ )z  , 

2 x 3x3y  x * 3x3z  x z * 


(4-2) 


3x 


where  * is  defined  as  follows: 
: 


♦z(x|x’) 


2 3z  (IA  + V ’ 


(4-3) 


where  I and  I„  are  given  by  (2-13)  and  (2-14).  The  electric  field  can 
A D 

then  be  calculated  by  using  proper  Hertz  potential  in  Chapter  2. 

For  a buried  horizontal  wire,  (Fig.  4-1),  the  incident  fields  have  to 
be  found  from  eq.  (4-1)  as  a function  of  horizontal  distance.  The  operator 
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Fig.  4-1.  Geometry  of  the  thin  wire  antenna  and 

scattering  problem.  All  distances  are  in 
meter,  where  = (antenna  length)/2, 

L2  - (wire  length) /2. 
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(4-2)  lias  to  lie  evaluated  by  finite  difference  methods.  In  order  to  avoid 

the  unnecessary  repel  I l ioii  oi  calculation  ol  ii^(x-x')  and  11  (x-x'i  lor 

different  x's,  we  shall  use  the  numerical  method  as  follows.  Instead  of 

computing  it  and  tt^  at  every  pair  of  specific  points,  we  shall  compute  in 

a subroutine  the  value  of  tr  (m*Ax)  and  tt  (m*Ax),  where  m ■ 1,  ....  M;  Ax 

X z 

is  a unit  distance  such  as  (antenna  length/number  of  divisions;;  M is  an 
integer  such  that  (M*Ax)  will  cover  the  maximum  distance  of  (x-x’>.  When 
the  values  tt^x-x')  at  points  other  than  mAx  are  required  interpolation 
methods  can  be  applied  easily  if  Ax  is  sufficiently  small. 

Fig.  (4-2)  shows  the  x-component  electric  fields  as  a function  of  kz 
for  fixed  kr  ■ 0.0833tt.  It  is  seen  from  the  phase  plot  that  the  electric 
field  behaves  like  a plane  wave  in  the  vertical  direction  as  we  would 
expect.  From  figs.  (4-3)  to  (4-7),  |Ex|  and  accompanied  phases  are 
plotted  as  a function  of  kr  for  various  fixed  values  of  kz.  These  magni- 
tudes are  reproduced  in  a three  dimensional  plot  in  fig.  (4-8) . 

With  the  incident  fields  so  determined  by  (4-1) , we  are  ready  to 
find  the  current  distribution  induced  on  the  buried  wire  (fig.  (4-1)). 

That  is,  instead  of  using  Vo*6(x)  for  Eox(£)  i-n  equation  (3-14),  the 
electric  field  found  in  (4-1)  is  used.  In  eq.  (3-14)  the  Hertz  potentials 

in  section  C of  Chapter  2 are  to  be  used  and  k , 1/e  art  to  replace  k,  e. 

t 

The  current  distributions  on  the  wire  for  different  locations  and  lengths 
are  shown  in  fig.  (4-9),  (4-10)  and  (4-11). 

Once  the  current  are  determined,  eq.  (4-1)  is  appiied  again  to  find 
the  scattering  fields.  For  a field  point  in  the  air,  i = 2,  j = 1. 

Results  are  shown  in  fig.  (4-12)  through  (4-16). 

As  a primary  check  of  the  algorithm,  we  check  the  results  with 
riciprocity  theorem.  That  theorem  states: 
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Fig.  4-10.  Current  induced  on  a horizontal  wire  buried  in  the  ground, 
a ■ (w/k)  , notation  is  defined  in  fig.  (4-1). 


(4-4) 


• J2dv'  - 


• J1dv'; 


where  indeces  are  used  the  same  manner  as  those  in  eq.  (4-2).  For  a 

particular  case;  f ■ .955  GHz,  e ■ 15  - 11.52,  L,  « .25X  , L_  ■ .129A  , 

X 0 / o 

r *=  .1751  , z0»  - . 3A  , and  h * . 33A  , (all  the  quantities  are  to  bo 
o/o  o 

refered  to  Fig.  (4-1)  and  Aq  is  the  free  space  wavelength).  The  integral 
on  the  right  hand  side  yields  a value  of  (-1.67353xlO~^-j3. 7149xlO-2), 
in  MKS  system,  and  that  of  the  left  hand  side  (-1.6727x10  ^-J3. 7056*10~2) . 


-65- 


References 


1.  Mei,  K.K.:  "Uni -moment  Method  of  Solving  Antenna  and  Scattering 

Problems",  to  be  published. 

2.  Stratton,  J.A.:  "Electromagnetic  Theory",  McGraw-Hill,  N.Y.,  (1941). 

3.  Harrington,  R.F.:  "Field  Computation  by  Moment  Methods",  The  Macmillan 

Company,  N.Y.,  (1968). 

4.  Sommer f eld,  A:  "Partial  Differential  Equation  in  Physics",,  Academic 

Press,  N.Y.,  (1964). 

5.  This  part  has  been  presented  at  1972  USNC-URS1  Spring  Meeting, 
Washington,  D.C.  April  13-15,  1972. 

6.  Somme rf eld.  A:  "Uber  die  Ausbreitung  der  Wellen  in  der  Drahtlosen 

Telegraphie",  Ann.  Physik,  81,  S.1138,  (1926). 

7.  Van  der  Pol.  B.:  "Theory  of  thr  Reflection  of  the  Light  from  a Point 

Source  by  a finitely  Conducting  Flat  Mirror,  with  an  Application  to 
Radio telegraphy,"  Physica,  _2  843  (1935). 

8.  Feynberg,  Y.L.:  "The  Propagation  of  Radio  Waves  along  the  Surface  of 

the  Earth",  English  Trans,  by  Wright-Patterson  Air  Force  Base,  Ohio, 
(1967). 

9.  Kuo,  W.C.  and  K.K.  ilei:  "Horizontal  Dipole  Arrays  over  Lossy  Ground", 

U.S.  AMERDC,  Fort  Belvoir,  Virginia,  Prepared  by  Electronics  Research 
Lab.,  U.  of  Calif.,  Berkeley,  Dec.  1971. 

10.  This  part  has  been  presented  at  International  IEEE  1G-AP  Symposium 
and  USNC/URSI  Meeting,  Boulder,  Colo.,  21-24  Aug.  1973. 

11.  Mei,  K.K. : "On  the  Integral  Equation  of  Thin  Wire  Antennas",  IEEE 

Trans.  AP-13,  373,  (1965). 


-66- 


Appendix  A.  Sommer f eld's  Integrals  for  an  Horizontal  Dipole 

For  a unit  infinitesimal  dipole  in  free  space  the  Hertz  potential  is 

given  by 


it 


1 -jfcR 
R 6 


(A-l) 


The  spherical  wave  representation  of  (A-l)  can  be  transformed  to  cylindrical 
wave  representation  by  Bessel's  transformation  to  obtain 


it 


e-j^R 

R 


f J (»r)e-vl*'h|  id*. 

Jo  V 

0 


(A-2) 


Where  all  distances  are  in  meter  and  shown  in  fig.  A-l;  J (*r)  is  the 

o 

Bessel's  function  of  the  zeroth  order,  and  v is  given  by 

v « A2-k2  (A-3) 

with  positive  real  part  to  insure  the  convergence  of  the  integral  and  its 
vanishing  at  |z| 

Now  Instead  of  considering  the  homogeneous  infinite  space,  the  haif- 
space  of  z < 0 is  the  earth  with  refractive  index 

er  * e'-j  e"/e°  , (A-4) 

where  c*  and  c"  are  the  real  and  imaginary  part  of  the  complex  permitivity; 
e"  is  the  loss  factor  for  the  ground,  which  is  dominated  by  conduction 
loss  in  high  frequencies.  Since  our  major  concern  is  application  in 
high  frequencies  hence  c^.  from  now  on  will  be  understood  by 

r » c * (c  '-jo/oj)  /e 
r J o 
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where  a is  the  conductivity.  The  nagnetic  permeability  in  the  ground  is 
taken  as  the  same  as  that  o 2 free  space  and  is  noted  as  Without  loss 

of  generality,  the  orientation  of  dipole  axis  is  chosen  along  the  x-direc- 
tlon  in  rectangular  coordinate.  The  Hertz  potential  it  is  related  to 
electric  and  magnetic  field  by  the  following  relations: 


E - V(V-n)  + k 7T  , 


(A-5) 


H ■ jui  lq(Vxtt)  for  z > 0, 


(A-6) 


where 


n « x n + z it  , 
x z 


- V(V-n£)  + k nE  , 


(A-51) 


H£  - ju)  e (Vxn£)  for  z < 0 . 


(A-6') 


Continuity  of  tangential  electric  and  magnetic  field  at  the  plane 


z*0  requires 


3 tt  r 3 it  _ 
x xE 

"x  b nxE*  9z  Sz 


(A- 7) 


3tt  9n  3tt  9n 

z zE  xE x 

% C ^zE’  9z  9z  9x  9x 


(A-8) 


Now  for  a point  source  above  the  ground,  we  may  imagine  that  there  is 
an  image  source  jn  the  ground,  which  produces  a similar  Hertz  potential 
as  in  (A-2) . The  solution  of  the  problem  for  z > 0 is  then  the  composi- 
tion of  contribution  from  the  primary  source  in  the  absence  of  the  ground 
and  secondary  source  induced  by  the  ground.  For  z < 0,  the  Hertz  potential 
will  be  represented  in  a similar  manner,  only  then  there  is  no  primary 
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source.  We  shall  consider  the  x-component  first. 


prim  . sec. 

X ■ ff  + Tt 

XX  X 


(A-S) 


OD  00 

- j Jo(Ar)e‘V!Z~hl  ^ dA  + j F^AjJ^Ar)* 


-v(z+h)d>  , 


OD 

"xE  ' / 


v_z-rh 

F2(X)Jo(Ar)e  dA 


(A-10) 


where  v£  ■ with  positive  real  part  and  F^(A),  F2(A)  are  weighting 

factors  yet  to  be  determined.  Using  (A— 9)  and  (A— 10)  in  (A-12)  we  find 
that 


Fi<*>  -7<'1  + 7T^>-  F2(i)  ' T T+^  ‘ 


And  the  solutions  are  then, 


-jkR  -jkR' 


x R 


OK 

/ 


i — +2  I J (Ar)e 


-v(;  f-h)  AdA 
v + v. 


(A-ll) 


2 ( , . VEZ-Vh  XdX 

"*E  ■ 7 J V*r)e  ttt. 


(A-12) 


where 


R2  = r2  + (z-h)2,  R’2  * r2  + (z+h)2. 


To  find  n and  ti  , note  first  that  there  is  no  primary  component 
z zb* 

for  tt  , since  the  dipole  is  x-oriented.  Note  also  that  there  is  ( 9 tt  /9x) 
Z X 

In  boundary  condition  (A-8) , and 

9t?  9 TT  9 TT 

x x 9r  , x 

9x  9r  9x  COS  Sr  ’ 
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? 


where  4 is  the  angle  between  x-axls  and  r,  hence  * must  contain  the 
factor  cos  4 and  will  be  constructed  by  J^(Ar)  instead  of  JQ(Ar).  Thus 


m 

1 ■ cos  4 / J, 

‘ J 1 


w,e"  cos  ♦ 


00 

A 


(Ar)e"v(l+h)«1(A)dA 


v_r-vh 

(Ar)e  6 *2(A)dA 


(A-13) 


(A-14) 


where  4^ (A)  and  4>2(A)  like  F^(A)  and  F2(A)  are  factors  to  be  determined. 
Solving  (A-8) , (A-13)  and  (A-14),  we  find  that 


2A2  V ‘ VE 

*1  — ;t  c we  •nJ  ’2 ' VE 


Therefore 


- 7-  cos 
k 


/ 


J^(Ar)e 


-v(z+h) 


v - v 


zE 


'~2  cos  4> 


CO 


(Ar)e 


VgZ-vh 


t v + v 


v - v_ 


r v + v. 


E A2dA 


A2dA 


(A~15) 


(A-16) 


Equation  (A-ll) , (A-12) , (A-15) , (A-16)  form  the  complete  set  of 
Hertz  potential  for  a horizontal  dipole  above  the  ground.  The  dependence 
of  both  source  and  observation  point  appears  only  on  the  exponential  terms 
of  all  integrals  in  these  equations.  These  exponential  terms,  however, 
are  crucial  to  the  convergence  of  integrals. 

Since  we  are  solving  boundary  value  problem  without  using  the 
knowledge  of  specific  property  of  parameters.  Hence  if  the  dipole  is 
located  in  the  ground  we  expect  to  obtain  a set  of  solutions  similar 
to  those  of  equation  (A-ll),  (A-12),  (A-15)  and  (A-16).  Only  now  k^,  k, 
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1/f  have  to  replace  k,  k^  v,  v£,  c respectively. 


Appendix  B.  Transformation  of  So— erf  eld 'a  Integral 

W«  wish  to  transform  sq.  (2-13)  to  ths  fern  of  eq.  (2-15). 


f 


JQ(Ar)e 


-v(z+h) 


v + Jk>/Tj 


dA 

v 


-Jk/c?  * 

Multiplying  both  sides  by  e we  get 

-jk/7^  r r -(v+jk v/^)z 

e 1 I - I JQ(Ar)e  1 e h XdX 


(v  + jk/^)> 


(B-l) 


Differentiating  both  sides  with  respect  to  z 

3 -Jk/c?  * f v -jk/c°  z 

~ (e  1 I)  - - | J (Ar)e“(t+h)  -*■  e 1 , 

d Z I O V 


-jut'  _:)k 

- e J e 


-ik/7j  z 1 "JkRt  *Jk/^  C 
• • I - - e I e dt 


R£ 


z 

•/ 


-Jk(Rj  + ((-*)) 


di  * 

r; 


(B-2) 


where  R'  and  R^  are  defined  the  sane  as  those  in  Chapter  2.  c is  the 
integral  limit  such  that  the  integral  converges.  It  is  proved  that  c 
may  be  chosen  as  -j  »,  since  at  this  limit  the  integrand  vanishes. 


{The  formula  and  figure  number  referred  to 
in  this  part  all  belong  to  this  part.] 


I.  Introduction 


Consider  the  electromagnetic  scattering  properties  of  dielectric 
cylinders  with  finite  cross-section.  If  the  dimension  of  the  scatterer 
is  of  several  wavelengths  or  larger,  the  problem  can  be  solved  by  geometrical 
optics  approximation.  On  the  other  hand,  the  well-known  low-frequency 
Rayleigh  approximation  is  suitable  to  solve  the  small  problems  with 
dimension  less  than  fractions  of  a wavelength.  For  the  problems  around 
the  size  of  a wavelength,  integral  equations  played  an  important  role  in 
computing  the  scattering  fields  [2],  [3],  [5],  (6).  Another  way  to 
calculate  the  fields  is  the  analytic  continuation  method  [4].  Recently 
the  blooming  techniques  of  the  Finite  Element  Methods  in  solving  boundary 
value  problems  enable  us  to  develop  a numerical  technique  to  tackle  the 
problem  with  great  versatility.  In  this  report,  the  Uni-moment 
Method  [1]  is  applied  to  calculate  the  scattering  fields  of  dielectric 
cylinders  of  arbitrary  cross-section  or  of  inhomogeneously  loaded  material. 

The  basic  technique  of  the  method  is  the  use  of  the  Finite  Element  Methods  in- 
side a mathematical  circle,  which  encloses  the  inhomogeneous  body.  The  fields 
outside  are  expanded  in  the  usual  cylindrical  expansion.  The  interior 
and  exterior  . roblems  are  first  solved  separately  and  later  coupled  by  the 
continuity  condition  of  the  E and  H fields.  The  versatility  is  increased 
greatly  by  introducing  the  method  of  "Inhomogeneous  Element".  The  advantage 
of  the  proposed  method  is  the  simplicity  and  efficiency  in  programming. 

Various  prob'ems  can  be  solved  by  the  same  program.  And  the  input  control 
for  the  shape  and  material  can  be  made  very  simple. 

The  following  discussions  will  include  both  E-wave  (E  i 0,  H =0) 

z z 

and  H-wave  (E  = 0,  H + 0)  scattering  of  dielectric  cylinders  of  arbitrary 
z z 

crossection,  or  of  inhomogeneous  material  such  as  l.unohorg  lens.  With- 
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out  any  further  effort,  scattering  of  multiple  cylinders  is  also  computable 
by  the  same  program.  The  incident  wave  is  remained  unrestricted  in  the 
method;  it  may  be  a plane  wave,  waves  generated  by  line  sources,  or  waves 
generated  by  other  kind  of  sources.  All  the  fields  computed  are  assumed 
to  have  the  conventional  exp(jwt)  harmonic  time  dependence. 


II.  The  Uni -moment  Method 


Consider  the  incident  waves  propagate  in  the  direction  normal  to 
the  z-axis — the  axis  of  the  dielectric  cylinder.  The  whole  system  is 
uniform  n the  axial  direction.  From  the  differential  form  of  Maxwell's 
equations,  we  know  that  E-wave  and  H-wave  are  separable.  Thus,  we  may 
investigate  these  two  kinds  of  fields  separately.  A circle  with  radius 
r = a is  drawn  to  enclose  the  entire  scatterer.  (Fig.  1). 


Fig.  1.  Construction  of  the  Mathematical  Circle 


inc  inc 

Let  the  incident  wave  be  represented  by  E for  E-wave  and  H 

z z 

for  4-wave.  Let  (or  H*)  be  tne  scattered  field  in  region  I (r  > a) ; 

£T  IX 

E * 'or  H ) be  the  field  in  region  II  (r  < a)  „ At  the  artificial  circular 
z z 

boundary  (r  = a) , the  tangential  components  of  E and  H fields  are  con- 
tinuous, i.e.,  for  E-wave: 


1 

i 

i 

i 


E l1  (a, A)  - E^n°  (a,4>)  + E*  (a,*); 


0 < <+>  < 2ir 


(1) 


1 


and 


3E. 

¥T 


ii 


(r  ,<!>) 


9E, 

r=a  " 3r~ 


inc 


r-a  + 17 


; 0 < J>  < 2tt 
r=a  — — 


(2) 
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4 

J 


"»'•»*  * 


•••*»•-  '*‘H**fc-  -*ver  **  --  *-* 


r*«-H  |H»1>  -- 


Or  for  H-wave: 

H^U,*)  - H*nc(a,$)  + H*(a.*);  0 £ $ < 2tt 

Z 2 Z 


d’) 


and 


3H 

i 

3r 


II 


,inc 


<r , 4>) 


r-a 


3r 


(r,<p) 


3H 

+ ~ 
r-a  3r 


(r,  «|> 


; 0 < <f>  < 2ir 
r-a  — — 


(2') 


Since  the  normal  derivative  value  of  at  a dielectric  boundary  is 

discontinuous,  it  is  convenient  to  have  the  artificial  circle  a little 

larger  to  make  sure  that  both  eqs.  (1')  and  (2*)  hold. 

E*(r,<f>)  [or  H*(r,<|))]  is  an  outgoing  wave  in  free  space.  The  general 
z z 

solution  is  expressed  by 

00 

E*(r,<J>)  [or  H1  ( r , d>)  ] = £ (k  r)  (A  cos  n <p  + B sin  n<f>)  (3) 

z z n non  n 

n=0 

(2) 

where  H is  the  Hankel's  function  of  second  kind  of  order  n, 
n 


k = W \/u  £ 

o o o 


A , B are  arbitrary  constants  to  be  determined, 
n n 

There  are  two  different  kinds  of  the  Uni-moment  Methods  which  work 
equally  well:  the  Dirichlet  type  and  the  Neumann  type.  We  shall  discuss 

these  two  methods  in  detail  separately. 


(A)  The  Dirich] et  type: 

The  main  idea  of  the  method  is  to  replace  E11  (a , 4>)  [or  (a, 4>)  ] by 

z z 

a linear  combination  of  a set  of  basis  functions  e^  ( 4> ) » n = 1»  2,  3,  ... 
of  the  continuous  functional  space  C[0,2tt].  That  is,  assume 
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(A) 


E1I(a,<fO  [or  H^1  (a,  <f>)  ] - £ cnen^ 
z n*l 


where  ( cf> ) are  assigned  known  functions  (e.g.,  polynomials,  or  sin  n<f» 

and  cos  n$) ; are  called  moments,  which  are  unknown.  Then,  by  the 

linearity  of  the  material  (i.e.,  c,  y are  independent  of  the  field  strength,) 
11  II 

E (r,<J>)  [ot  ( r , 4>)  ] can  be  expressed  by 

z z 


E^Cr,^)  (or  H^1  (r , <J>)  ] * £ cn*n(r,+)  . (5) 

n=l 

In  eq.  (5),  ^(r.if))  is  the  solution  of  the  following  system: 


E-wave: 

V2^(r,4>)  + *V  r (r , <l>) ^n<r,4>)  = 0 


(6) 


H-wave : 

v-t 


—77  r<-  (r , 4>)  ] + (r,4>) 

e I r , <p)  n on 


= 0 


subject  to  the  boundary  condition 


~ en(<f>)  ; 0 <_  <J>  <_  2tt 


(7) 


(8) 


where 


(-,.(r,<|>)  = e(r,4 >)/c 
r o 


The  numerical  solution  of  ij;^(r,<f>)  are  generated 

Methods  which  will  be  described  in  the  next  section. 

ip  (r,<p)  obtained,  it  is  possible  to  calculate  — — (r 
11  i)r 


by  the  Finite  Element 
From  the  values  of 


♦) 


r -a 


Hence 
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)E  3H  * dVn 

iT"  <'•+>  r-a  lor  ST-  r-aJ  ’ £ cn  ~ <*•♦> 


r*a  . 


From  eqs.  <3 ),  (4),  and  (9),  plugg  into  eqs.  (1)  ind  (2)  [or  eqa.  (1', 
and  (2')]  we  get 

QO  00 

5]  c e (*)  • EinC(a,<f>)  [or  H*nc(a,4>)]  + £h!;2)  <kQa>  (An  cos  n<J>  + Bn  8in  n*J 
n«l  n n 2 2 n-0  n ° 


• dfc 

£ cn  s2  (r>*>  * jt-  <'•♦>  ‘°r  <r-*>  1 

n=l  r«=a  r=a  r-a 


/ O'i  1 

+ Y*  k Kv  ; (k  a)  (A  cos  n<j>  + B sin  n<f>) 
o n on  n 

n=0 


The  unkncwn  coefficients  An's,  B^s,  and  Cn's  can  be  solved  from  eq.  (10) 
and  eq.  (11).  Approximate  solutions  are  obtained  by  either  point-matching 
or  weighted  least-square  method.  It  is  more  convenient  if  we  choose  en(<fr) 
to  be  trigonometric  functions  sin  nif>  and  cos  n$.  If  this  is  the  case, 

00 

EII(a,4i)  [or  HU(a,$)]  = £ (C  cos  n*  + Dn  sin  n<t>)  . (12) 

z 2 n=0 

Let  (r,<p)  be  the  solution  of  eq.  (6)  [or  eq.  (7)]  with  the 

boundary  condition 


^(a,<J>)  = cos  n<|>  (13)  , 


i|>  (a,$)  = sin  n<p  . 
n 


Similar  to  eq.  (10),  we  have 
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(15) 


E<c  cos  n<J>  + D sin  n<p)  « Einc(a,$)[or  Hinc(a,<f>)] 


n=0 


+ E H^(k  a)  (A  cos  n<j>  + B sin  n<[>)  . 
n*0  on  n 


Similar  to  eq . (11),  we  have 


<*>  3ip 

£ <c„  IT  <*•♦> 

n=0 


34) 

+ D r— ^ - (i;4») 
r=a  n 3r 


3E 


.inc 


J = 


r-a  3r 


(r*<t>) 


3H 


inc 


for 


r-a  — ar'~  (r-*)|r-a1 


oo 

r~>  ( 9 ) * 

+ E k " (k  a)  (A  cos  np  + B sin  np)  . 


n=0 


on  o n 


(16) 


So,  from  eq . (15) , 


C - A H(2)(k  a)  + fC 
n n n o n » 


(17) 


D = B H(2)(k  a)  + fS 
n n n o n * 


where 


2tt 

fC  = - f 

n 7T  I 

•fn 


Einc(a,4i)  [or  H ^nc(a,ij>)]  cos  n<J>  d<j>  , 
z z 


(18) 


(19) 


2n 

fs  . a r 
n ' i 


E^nc (a , <J>)  [or  H*nc(a,<ji)]  sin  n<f>  dij> 


(20) 


From  eqs.  (16),  (17),  and  (18), 


£ >An(H!S2)  <ko“>  7T 

n=0 


dp 


-a  - ko^2>  ' <koa>  cos  n*>  + BX"<koa)  TT  <r-*> 


(2), 


dp 


t ^e^110 

“ koHn2)  (koa)  S±U  "■  JT (*.♦) 

°°  3^C 

- E Iff  (r,«fr) 


3H 

r=a  ^°r  3r 


inc 


(r,<j>) 


r-a 


n=0 


n 3r 


dp8 

+ fS  ~ (r ,p) 
r=a  n 3r 


r=a> 


(21) 


r-a 
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The  above  summation  is  truncated  to  N terms.  Following  the  analysis  given 

by  Wilton  & Mittra  (4)  and  the  numerical  data  available,  N is  chosen  to  be 

a number  slightly  larger  than  kQa.  Then,  multiplying  both  sides  of  eq. 

(21)  by  sin  n<(>,  (n  ■ 1,  2,  . . . N)  and  cos  n<f>  (n  ■ 0,  1,  2,  ....  N)  an' 

integrating  from  0 to  2tt,  we  end  up  with  a set  of  linear  equations  which 

can  be  used  to  solve  for  A ’s  and  B ’ s.  Note  that  only  the  right  hand 

n n 

side  vector  is  affected  by  the  incident  fields. 


(B)  The  Neumann  type: 

The  idea  is  quite  similar  to  the  Dirichlat  type.  The  only  difference 
is  to  assume 


3E 


i.1 


3r 


9H 

r-a  Ior  5T 


II 


(r  ,4)) 


) «=  53  C e (4>) 
r*a  ~ n n'r/ 
n=l 


(21) 


where  C , e are  defined  as  before.  Following  the  same  argument  as  that 
n n 

in  the  DIrichlet  type,  eq.  (5)  holds,  and  i^(r,4>;  is  the  solution  of  eq.  (6) 
or  eq.  (7)  subject  tc  the  boundary  condition 


r 

3r 


(r,<J>) 


_ „ = e ($) ; 0 <_  (j»  _<  2tt  . 

i-a  n — — 


(22) 


The  numerical  solution  of  1//  (r,4>)  are  generated  by  the  Finite  Element 

n 


methods.  Hence, 


E^1  (a,<J>)  [or  H^a,*)] 
z z 


E C ip  (a,4>)  • 
. n n 
n=l 


(23) 


From  eq.  (21)  and  eq.  (23),  similar  results  as  eq.  (10)  and  eq.  (11) 
are 


n=i 


CJn(a,<P)  = E*nC(a,<mor  H*nC(a,$) 


E H 


(2) 


n-0 


(k  a) (A  cos  n$  + B sin  n4) 

o n n 

(24) 
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and 


® 3E 

L C e ($)  - r-5 
, n nVT'  3r 
n*l 


inc 


(r,4>) 


dH 

r-»  lor  IT 


Inc 


(r  ,$) 


r-a> 


+ k H ^ (k  a)  (A  cos  + B sin  n«#>) 
*— f o n o n n 

n-0 


(25) 


The  unknown  coefficients  A 's,  B ' s,  and  C 's  are  then  solved  from  eqs. 

n n n 

(24)  and  (25)  in  a similar  way  as  that  of  the  Dirichlet  type. 

Several  advantages  have  been  found  in  the  methods  discussed  aboveJ 

(1)  The  artificial  circle  enables  us  to  solve  the  interior  and  exterior 
problems  separately.  Only  the  interior  problem  involves  inhomogeneous 
materials  of  various  scatterers. 

(2)  In  case  the  boundary  conditions  of  the  interior  problem  is  expanded 
into  Fourier  series,  the  number  of  terms  in  the  infinite  summation  to  be 
truncated  can  be  determined  a priori. 

(3)  Since  the  interior  and  exterior  problems  are  coupled  at  the  circle, 
the  Hankel's  function  and  its  derivatives  are  computed  only  once  at  r ■ a, 
the  radius  of  the  artificial  circle. 


(4)  Host  of  the  computations  do  not  involve  the  incident  fields.  Calcu- 
lation for  various  Incident  waves  add  only  little  effort. 

(5)  The  versatility  is  developed  in  the  Finite  Element  Methods  which 
can  solve  the  irregular  problems  relatively  easily.  The  method  of 
"Inhomogeneous  Element"  enables  us  to  compute  various  problems  by  a 
single  program  and  to  simplify  the  input  control. 
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III.  The  Finite  Element  Methods 


It  is  impossible,  for  most  cases,  to  solve  4/  (r,$)  analytically  for 

n 

differential  equation  eq.  (6)  or  eq.  (7)  with  boundary  conditions  eq.  (8) 
or  eq.  (22).  The  Finite  Element  Methods  for  solving  partial  differential 
equations  has  been  investigated  by  many  civil  engineers  and  some  authors 
in  Electrical  Engineering  such  as  Silvester  and  Wexler  [/],  [8].  It  can 
be  shown  that  the  solution  of  eq.  (6)  subject  to  the  Dirichlet  boundary 
condition  eq.  (8)  is  equivalent  to  solving  the  function  which  gives 
the  stationary  value  of  the  variational  formula  I, 


I Vi|/  |2  - k2e  tj>^)ds 
1 n1  o r n 


(26) 


The  integration  region  is  the  entire  disc  enclosed  by  the  artificial  circle 
(r  1 a)  • 

With  some  manipulations,  the  variational  formula  for  eq.  (7)  subject 
to  the  Dirichlet  boundary  condition  eq.  (8)  Is 


2 2 

- k j )ds 
on  • 


(27) 


The  region  S of  the  integration  is  the  same  as  that  of  eq.  (26). 

Following  the  procedures  suggested  by  Wexler  [8]  , th>?  variational 
formula  for  the  Neumann  boundary  condition  eq.  (22)  can  be  derived  from 
eq.  (26)  and  eq.  (27).  That  is,  the  formula  for  eq.  (6)  subject  to  the 
boundary  condition  eq.  (22)  is 


I 


( I Vip  I2  - k2e  ij>2)ds  - 
1 n 1 o r n 


^nb($)d£ 


(28) 
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where 


b($) 


(r,4>) 


r*a  * 


(29) 


The  region  of  integration  S is  the  entire  disc  enclosed  by  the  artificial 
circle  (r  a)  ana  C is  the  contour  of  the  circle.  The  direction  of  the 
contour  integral  is  counter-clockwi  ,e. 

Similarly,  the  variational  formula  for  eq . ,7)  subject  to  the  boundary 
condition  eq.  (22)  is 


I 


2 2 

- k>  )ds 
o n 


(30) 


where  b(<J>)  is  defined  in  eq.  (29).  The  region  and  contour  of  integration 
are  the  same  as  that  of  eq.  (28).  The  integrand  of  the  contour  integral 
in  eq.  (30)  should  be  divided  by  e^.  But  since  the  artificial  circle  is 
larger  than  the  scattcrer,  is  unity  at  the  contour. 

The  region  S is  then  divided  into  tn3ny  small  triangular  elements. 

The  way  to  construct  the  triangles  is  very  flexible.  The  only  restrictions 

are  convergence  and  prograninability.  It  is  not  necessary  to  consider  the 
inhomogeneity  of  the  material  in  constructing  the  elements.  The  method 

of  "Inhomogeneous  Element"  will  handle  the  inhomogeneity  or  irregular 
boundary  of  the  scatterer.  Fig.  2.  is  an  example. 


(A)  Dirichlet  problems: 

Assign  a number  to  each  node  of  the  triangles  in  the  circle.  Let 
ip  be  the  value  of  the  function  ip  at  the  i-th  node.  Let  i = K+l,  K+2, 

K-t-K'  be  the  boundary  nodes  on  the  artificial  circle,  and  i • 1, 

2,  3,  K be  the  nodes  in  the  interior  of  the  circle.  Approximate 

the  function  in  each  triangular  element  by  a linear  function  which  passes 


g.  2.  Example  of  Triangular  Elements 


symmetric 

axis 


°VV 


3.  Inhomogeneous  Element 


through  the  functional  values  at  three  nodes.  Then  integrate  eq.  (26)  or 
eq.  (27)  through  that  elenent.  Add  the  contributions  of  all  elements 
together  to  get  the  global  integration  of  I,  which  is  a quadratic  function 

of  * t i - 1,  2,  ...  K,  «C+1,  K+2 K+K\  To  make  the  functional 

1 stationary,  take  derivative  of  1 with  respect  to  ^ and  let  it  be  zero, 
i.e.,  ~ * 0,  i - 1,  2,  ....  K.  Since  I is  qw  *tic  function  of 


the  above  equations  become  a system  of  linear  equations,  which  can  be 
written  in  the  matrix  form: 


(31) 


The  matrices  S and  B are  obtained  from  the  integration  ot  all  elements.  The 


vector  [<1^,  1 1>2> 


T 

1 is  to  be  solved.  The  vector  t^K+1> 


n+2  ' 


} ijj  ,]T  is  an  array  of  the  boundary  values  which  are  known, 


The  unknown  vector  can  be  solved  by  Gaussian  elimination  of  eq.  (31), 
or  the  block  by  block  method,  or  a sparse  matrix  algorithm.  Once 

9ili 

iji  , i = 1,  2,  3,  .....  K are  known,  the  normal  derivative 


, which 


is  required  by  section  II,  is  readily  calculated. 

The  rest  of  the  problem  is  to  construct  the  matrix  S and  B.  This 
requires  the  integration  in  every  element.  If  the  material  is  uniform  in 
an  element,  the  formula  of  the  element  has  been  derived  by  Silvester  (7). 
We  shall  discuss  the  cane  when  the  dielectric  constant  in  the  element 
is  a continuously  varying  function,  and  the  case  when  the  element  contains 


j 
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two  different  materials.  This  technique  is  called  the  "Inhomogenecus 
Element"  method. 


(1)  The  triangle  contains  two  materials: 

For  this  case,  we  approximate  the  dielectric  boundary  in  the  element 
by  a straight  line  as  in  Fig.  3.  This  approximation  is  valid  whenever  the 
scatterer  surface  is  quite  smooth  or  the  element  is  quite  small. 


(i)  E-wave: 

Since  E and  its  derivatives  are  continuous  accross  a dielectric 
z 

boundary,  a linear  function  is  used  as  a trial  function. 


4/  * ax  + by  + c 


02: 


(33) 


(34) 


(35) 


(36) 
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*j  " a2XJ  + b2yj  + C2  ) 

*k  " alXk  + bl\  + C1  y • 


(42) 


(b)  Continuous  at  both  sides  of  the  line  segment  between  point  1 and 
point  2,  i.e. , 


^alxl  + blyl  + cl^ 


(a^  + b1y2  + C;L) 


(a2xi  + b2y1  + c2) 
(a2x2  + b2y2  + c2) 


(43) 


(c)  The  normal  derivatives  at  the  line  segment  between  point  1 and  point 

I 

2 are  related  by  the  following  equation: 

1 ..  ,1  ^ I „ .2  ' 

— V <i>  • n = — ViJj  • n . 

t r 


[This  equation  is  derived  from  the  continuity  of  the  tangential  components 
of  the  electric  field.]  So, 


Ib1(x1-x2)  - a1(y1-y2)]  = ~ [b^-x^-a^-y^ ] 


Combine  eqs.  (42),  (43),  and  (45)  into  a matrix  form: 


N a = M \p 

where 


(45) 


(46) 
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and 


" H2  ! * 


(54) 


From  eq.  (27), 


, ■ // 


ijk  A 


ijk 


S 


ijk 


( 7-  \**\2  - kVjdxdy 


I I7*2!2  - kl*2 


(55) 


- ff  iv*2i2  - ty2nxd,  + /jf(  i_ 

2 412k  1 

jv/l 2- 

Fro*  eq.  (40)  and  eq.  (53), 

r~ 

1 .2  2 

k x 

e2  o 

.2 

- k xy 
o 7 

- k2x 
o 

~ n2i2-*2/ - i.,b,cj 

.2 

- k xy 
o 

1 2 
7~  - k y 
e2  ° 

-k2 

o 

- k2x 

L.  ° 

.2 
- V 

-k2 

o 

• 

1 ,2  2 

k x 

e2  0 

- k2xy 
0 

- k2x 
o 

T T 

- ♦ h\ 

.2 

- k xy 

1 ,2  2 

k y 

$2  o' 

•» 

- V 

- k2> 
0 

w 

• koy 

-k2 

o 

Similar  result  can  be  obtained  for  — I VP1 

E1 

I2  - k2*1*. 

]dxdy . 


( 

b 

c 


i 

J 


S11  * 


WTHWlltl  WJUllimil  lip,>  v^,.  9W* 


"i  111 


■ ' ■ \jk  ■ -I|R*iiHi + h*w* 


where 


e2  *poo"poo)"ko*p2<Tp2<P 


‘ ko(pxrpii> 


and 


Ty  * 


■ ko(p10-p10) 


Cj  poo  " kop20 


.2- 

- koplI 


- kop10 


e2  ^poo  roo)-ko(po2-po2) 


■ ko(pol-pol) 


koPll 


1 - ,2- 

— p - k p _ 

oo  o o2 


- k2?  . 

o ol 


- kop10 


,2- 

* k p , 
cv  ol 


- k2? 

o oo 


(56) 


ko(p10‘P10J 


ko(porpol> 


k (p  -p  ) 
o oo  *oo 


(57) 


(58) 


The  definitions  of  p and  py?  in  eqs.  (57)  ind  (58)  are 


// 


v C 


vc  x y dxdy;  v • 0,  1,  2,  ....  c - 0,  1,  2,  


VC 


ijk 

// 

AJ2k 


xVy^dxdy ; v - 0,  1,  2 C - 0,  1,  2,  


(59) 


(60) 


llic  integration  of  eq.  (59)  and  eq.  (60)  can  be  performed  in  several 
ways  (91. 

(2)  The  dielectric  constant  is  a continuously  varying  function: 
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Fig.  4. 


Let  eqs . (32)  through  (38)  all  hold  and  the  same  notations  are  used, 
(i)  E-wave: 

Approximate  er(x,y)  in  A^  with 
er(*.y)  - Cx  + ny  + c . 

Let  ev  “ er<xv’yv)  Whore  v - i,  j,  k . 

Since  is  a known  function,  C,  n,  C «re  solved  by  eq.  (63). 


r- 

e 

.,-1 

ei 

n 

■ N 

eJ 

ek 

L J 

From  eq.  (26) 


■ S 


ijk  A 


( | V*2 | 2-k2er^2) dxdy 


T -1  T -l 

*(N  V T3  (N  l)* 


(61) 

(62) 


(63) 


(64) 
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where 


f\J‘ko(tp30+',p21+cp20)  -ko<tpSl+np12',cpn)  ■koUp20t  npU+  tplO)" 

-ka^p21+np12+tpn)  poo-ko(tp12+,po3+tpo2)  'ko(Cpl]*npo2+Cpol) 


■ko('.P20t"PuHP10) 


-ko(tpll+np„2^pol) 


'ko<eplOVf,polHpoo) 

(45) 


PvC*  v * °»  lp  2»  3,  ....,  ; ■ 0,  1,  2,  3,  are  defined  in  eq.  (59). 


(ii)  H-wave: 

Approximate  In  t(Jk  with 

vbr  * tx  + ny  + c 


Let  the  assumption  of  eq.  (62)  hold, 


£ 

n 

■■  N**1 

i/Ej 

c 

1/Ek 

— — * 

_ J 

Then  from  eq,  (27) 


< |Vij/J2-k^2)d*dy 


T -1  T -1 
* (N  ) T4(N  A)* 


(66) 


(67) 


(68) 
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■ 1 ipj._ni.-jni 


where 


"<p10*"Pol  + '-ko»M 

' kopn 

1 

0 

1 

T4* 

- *1*11 

fpl0+"pcl+<‘kopo2 

■ Vol 

(69)i 

kop10 

-kV 

0 ol 

- kopoo_ 

p^,  v • 0,  1,  2,  ....  C « 0,  1,  3 ...  arc  defined  in  oq.  (59). 

(3)  Neumann  Problem : 

If  we  uiuat  the  boundary  condition  to  be 


5r 


Hi)  . 


(70) 


the  varlatlonel  fonaulae  eq.  (28)  or  eq.  (30)  should  bo  used. 

The  first  part  in  either  fonaulae,  i.e.,  the  suciace  integral,  is 
exactly  the  sane  aa  that  in  eq.  (26)  or  eq.  (27).  The  onJy  two 
differences  are: 


(a)  At  the  boundary  points  i-K+1,  K + K',  the  function  i is 

unknown,  but  we  restrict  to  have  the  following  property. 


31 


0 


for  i - 1,  2, 


K,  K + 1, 


K + K' 


(71) 


In  sitrlx  fors,  eq.  (71)  can  be  written  by 


~*1 

~ o ~ 

*2 

• 

m 

« 

• 

• 

• 

0 

• 

^K+K' 

* 

K cleaents 


(721 


K'  elements 
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(b)  The  vector  B In  eq.  (72)  ia  derived  f xm  the  line  Integral  of  ».ha 
variational  formula.  We  rhall  dlaeuaa  the  line  Integral  In  detail  because 
It  requires  aoae  technique*  In  aplina  analyela.  It  la  eaaentlal  to 
evaluate  the  line  integral 

[a,+)bU)d*  . (73) 


t*  - [ *0 


the  reault  of  eq.  (73)  la  in  tenia  of  a linear  function  of  4^,  1 - k+1 , 
K+K'.  Then 


a*i+K  1 


where  ia  the  i-th  element  of  B in  eq.  (72) 


’l-l 


■i+l 


Fig.  5. 


Let  *i+l  *i-l  " 6 


(7  <>) 


Let 


0 ■ $ - 


(75) 


The  Idea  is  to  approximate  tjt(a,t)  by  a quadratic  polynomial  of  h for  the 
interval  {-  j J 


tl»(a,$)  - ty(a,8)  ■ a02  + be  + c;  - 


(76) 
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And  tha  polynomial  aatiafiaa 


♦(•.-6)  - • ad  - b6  ♦ c 


<K«,0)  - * a*0  b<0  -f  c ■ c 


*(a,6)  ■ » *i S b«  f c 


Solw  uq . (77),  wc  have 


^i+K  ^ 


26 


^i-fl-HC  " *i-l+K 


26 


c - * 


i-HC 


Fro*  aq.  (76)  and  aq.  (78),  tr*  hava 


Lat'a  partition  aq.  (73)  by  latting 

A+ f 

li  m J ♦(a,6)b(6)dq  . 


♦i  ’ 2 


Changing  of  variable  6 to  8,  and  usa  aq.  (79), 


l‘-J> 


32 


a - 

6 


62  8 


(77) 


(7£) 


*(a,8)  - P141C  (1  - 2 * * *i-H*  * 7T  ” 26^  * *i+l+K  * 2 * 26  * * *79) 

0 26  26 


e2  . e 


i+K  '*  ,2  ' ’ ri-l-HC  * „;2  " 26  * * *i+l+K  ( „,2  * 26^b(V*)d* 


26* 


26* 


" Vi+K  * Vi-l+K  * Vi+l+K 


(80) 
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igtfBKt’it; 


1 


where 


6 

il 


(1  - ~ ) b<« >8)d0 
6 


(81) 


r*«-4- 

l i_  ** 


2J  )b(ii-H))de 


(82) 


i'fl  < $ ♦ jf 


(83) 


The  totel  integration  I ' becomes 


K'  K* 

X'  " ^ *i  " ^ <Ri*i-HC  * + 1-HC^ 


**i  * Si*1  + Tl"1^  *i-HC  ‘ 


(84) 


From  eq.  (84)  and  eq.  (73),  we  get 


VJ(Ri  + Sm  + IH)  • 


(85) 


It  is  important  tc  note  that  the  numerical  solutions  of  R^,  S^,  and 
froei  eq.  (81),  eq.  (82),  and  eq.  (83)  may  be  very  inaccvtate  when  6 
is  very  snail.  This  is  due  to  the  fact  that  the  integration  cosms  out 
with  -I  very  saiall  value  ia  divided  by  a email  vt»lue.  Round  off  error  is 
very  significant.  One  posalole  way  to  overcome  this  defect  is  to  expand 
K^+9)  in  Taylor  series  expansion  around  4^  before  integ  at  ion.  That 


is 
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""yjt-  2"~  “ ■-  ■ '•  .j.  * jjjgflyri y '*■•_»;.  ;j;.  r • '*4*T" 


IV.  NtaserlcaJ  Results: 


Bated  on  the  theories  ditcuaaad  in  this  report,  a FORTRAN  program 


has  bean  written  to  compute  the  expansion  coefficients  and  far  Melds  of 


the  scattered  waves  of  dielectric  cylinders  with  arbitrary  shape  or  in* 


homogeneity  Jnciuunted  by  a plane  wave.  A CDC  6400  computer  has  been 


used  to  run  the  prograa.  The  validity  of  the  program  is  demonstrated  by 


the  automatic  error  checking  for  dielectric  circular  cylinders.  For  this 


case,  the  prograa  prints  out  both  the  numerical  and  exact  values  of  thu 


Finite  Element  solutions  and  the  resulting  scattering  fields.  Also,  the 


results  for  dielectric  elliptical  cylinder  with  ka  ■ 1.2,  kb  ■ 1.0198, 


cf  • 2.0,  E-wave,  mat chew  perfectly  with  that  given  by  Wilton  t Nittra 


(4|.  Bee idea,  the  dielectric  ring  problem  is  coopered  with  Richmond's 


result.  As  ia  shown  in  r'ig.  6.,  the  solution  of  the  uni-ament  met  hoc  is 


closGr  to  the  true  solution  than  that  given  by  Richmond  1 3). 


Two  kinds  of  error  way  exist  to  affect  the  final  result  of  the 


scattering  fields:  the  truncation  error  of  the  Infinite  series  in  the 


Uni-moment  metheds,  and  the  error  due  to  the  Finite  Element  Methods. 


There  is  not  muen  thing  that  we  can  do  for  the  first  etror  since  the 


terms  to  be  truncated  are  determined  n priori  according  to  t!«e  size  of 


the  scstterer.  The  contribution  of  the  second  error  is  presented  in 


tables  1 through  3 . It  can  be  seen  from  these  tables  that  when  the 


normal  derivative  values  at  the  circle  are  relatively  small  for  some  modes. 


the  relative  error  of  the  finite  element  solution  may  be  very  large.  But 


this  affects  the  final  result  very  little  because  of  its  small  value.  An 


thu  normal  derivative  values  for  some  modes  being  large,  the  interior 


syrtem  reaches  resonance  condition.  The  error  for  this  mode  is  eo  large  as 


we  could  expect.  Fortunately  enough,  the  physical  systems  jre  mostly  irregular. 


(or  which  the  resonance  conditions  arc  usually  very  complicated.  The  con 
vergence  property  of  the  method  is  shown  in  Fig.  ?.  It  is  easily  seen 
thav  the  numerical  results  converge  to  the  true  solution  as  long  as  the 
else  of  the  elements  bscoming  smaller  and  smaller.  The  computation  time 
as  a function  of  the  number  of  unknowns  is  plotted  in  Fig.  6.  The  number 
of  unknown  is  inversely  proportional  to  the  square  of  the  element  site. 
Some  other  results  are  also  shown  in  Fig.  9 through  13.  in  which  several 
Interesting  geometries  were  investigated  to  show  the  capability  of  the 


t 


V.  Conclusion 

The  analysis  and  results  presented  In  this  part  are  solid  evidence 
that  the  "uni-moment  method"  is  a valid,  economical,  numerical  method  to 
compute  scattering  problems  involving  inhomogeneous  media.  We  are  conflden 
that  it  will  be  successfully  applied  to  many  interesting  problems  in  field 
theories . 
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Table  l.w 


h™  1 unit,  e ■ 2,  a ■ .I.\,  A n . 1125A  ka  * .89  E-wave 


a.  Error  of  FEM  Solution 


b.  Expansion  coefficient 


en 

error 

Const. 

.3 

2.8*10  3 

2.9*10“2 

cos  4 

1.3*10“-* 

4.4*10“2 

sin  4 

1.9*10”3 

4.4*10~2 

cos  24 

2.4*10~3 

l.OxlO-1 

sin  2$ 

2.2*10-3 

l.OxlO-1 

\> 


Computed 


-2 


-9.61x10  -j .295 


Exact 


-2 


-9.65*10  -j .295 


1 


3.09*10~2-j4.8*10“4 


3 . 10*10“2- j 4 . 8*10-4  } 


too  small 


c.  Far  Field  Pattern 


n 

4 

0° 

45° 

— 

90° 

135° 

180° 

lEJ 

Computed 

.281 

.289 

.310 

.331 

.340 

Exact 

.282 

.290 

.310 

.332 

.341 

*Ez 

Computed 

-110° 

-109° 

-108° 

-107° 

-107° 

Exact 

-110° 

-109° 

-108° 

-107° 

-107° 

(*)  In  tables  1 through  3,  the  notations  are  defined  as  f.  ; following: 


h • relative  element  si.:e.  1 unit  ; gg  A 
a » radius  of  the  dielectric  circular  cylinder. 
A ■ radius  of  the  artificial  circle. 

21i.  «,«  r it 

O'  < * < 360' 


■“  'SI 


«r°r  >1  E I §f  (A.tp  - ff  (A.^j)  | ]/[  £ | (A.^Xl 

where  prime  means  computed  value. 
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Table  2 


h ■ 2 units,  er  ■ 2,  a ■ .2A,  A »»  . 2250A,  ka  ■ 1.78,  E-wave 


a.  Error  of  FEM  Solution 


e 

n 

erxor 

const. 

7.4xlO-3 

coa  t 

Q3JHI 

sin  $ 

coa  2+ 

7. 3*10~3 

sin  2$ 

5.2xl0“3 

cos  3$ 

sin  3$ 

b.  Expansion  coefficient 

Computed  Exact 


Aq  -.403-j .491  -.406-J.491 


. 494- j. 130  ,496-j .132 


A2  2.19x10""4+J 2 .89*10  J 4.58xl0'V)3.03xl0“2 


too  3mall 


Exact 


0° 

45® 

90® 

135® 

180° 

.274 

.342 

.613 

.974 

1.15 

.275 

.342 

.614 

.978 

1.15 

-174.71® 

-155.5" 

-131.11® 

-120.55® 

-117.76" 

-175.03® 

-155.76" 

-131.32® 

-120.66" 

-117.84® 

Table  3. 


h « — units,  c * 2, 
3 r 


a.  Error  of  FEM  Solution 


■ 

error 

“x# 

conat. 

4.7*10"2 

1.8xlO~2 

cos$ 

3.39 

9.72 

ain+ 

3.35 

9.72 

cos  24 

2.2xl0“2 

5.6xlO’2 

sin  24 

1.88X10"2 

5.6xl0“2 

cos  34 

gggH 

4.4xl0“2 

ain  34 

1.7xlO~2 

4.4xl0“2 

a • .41*  A * .4331  ka  ■ 3.55,  E-w»ve 
b.  Expansion  coefficient 


Computed 


Exact 


Computed 


Exact 


Field 


0° 


.476 


.482 


-112* 


-114° 


Pattern 


45* 


.453 


.457 


-3.8° 


.5 


1 ■'  ■ * 1 

Computed 

Exact 

-.625-3 .484 

-.632-J.482 

.836-31. 55 

.828-j 1.56 

1.13+3 .992 

1.151-3 -989 

-.2131-3.023 

-.217+3.024 

*o 

°. 

1 

i 

135“ 

180“ 

1 

.706 

1.93 

4.18 

.713 

1.94 

4.21 

44.4° 

-152° 

• 

**> 

H 



43.5° 

-152“ 

-143° 

Echo  width/ X 


30  60  90  120  150  180 


<>  (Degrees) 


o*  ?.  X 

Com  A:  E-wom 


ComB 

£Jk  «,=  «/«.=  2 0 

ComC  E-wOM 
ComOH-mm 


tan*  rodiu«  * 25X 
Outor  roAul  *.30X 

<r  * </«•*  ^ 0 

C«t«  E 


Fig.  9.  Scattering  by  Various  Scatterers. 
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Echo  vidth/X 


6C  120  180  240  300  360 

« (Degrees) 

rig.  11.  Scattering  by  r^uare  Cylinder. 


Echo  width/ A 


(1)  R.  R.  Mel,  "Uni'MMAt  Method  of  Sol  vine  Antenna  and  Scattering 
Problems."  to  be  published. 

(2)  R*  R.  Mel  and  J.  ti.  Van  Sladel,  "Scattering  by  Perfectly •Conducting 
Rectangular  Cylinder a,"  TREK,  AP-11 , Kar-h  1961,  p.  185. 

(1)  J.  H.  Richmond,  "Scattering  by  a Dielectric  Cylinder  of  Arbitrary 
Croaa  Section  Shape,"  IEEE.  AP-11,  May  1965,  p.  144  and  AP-14, 

July  1966,  p.  460. 

(4]  D.  R.  Wilton  and  R.  Mlttra,  "A  Nev  Numerical  Approach  to  tha 
Calculation  of  Electromagnetic  Scattering  Properties  of  Two- 
Dlmaoeional  Jodies  of  Arbitrary  Croaa  Section,"  IEEE,  AP-20,  May  1972, 

p.  110. 

IS)  N.  Okamoto,  "Matrix  Formulation  of  Scattering  by  a Homogeneous 
Cyrotroplc  Cylinder,"  IEEE,  September  1970,  AP-18,  p.  642. 

(6)  T.  C.  Tong,  "Scattering  by  a Dielectric  Rectangular  Cylinder," 

1971,  0-AP  International  Sympoelum.  p.  79. 

(7)  P.  Silvester,  "High- Order  Polyncw'al  Triangular  Finite  Elements 
For  Potential  Problems, " Int.  J.  Eng.  Scl..  1969,  vol.  7.  p.  849. 

(3)  A.  Wexler,  "Computation  of  Electromagnetic  Fields,"  IEEE. 

MTT-17,  August  1969,  p.  416. 

(9)  Strang  and  Fix,  "An  Analysis  of  the  Finite  Element  Methods," 

Prentice  Hall,  1971. 


-40- 


